function result = plexIntShowAnalog(option)
% Plot Plexon analog data for the Snow Dots plexonIntegration demo.
%
%   showPlexonAnalog() uses the Plexon Matlab client SDK to access analog
%   online, during the Snow Dots plexonIntegration demo.  It expects be
%   called periodically by a  dotsComputableFunction object on the
%   server-side of dotsTheComputablesManager (i.e. on the Plexon machine).
%
%   The behavior of showPlexonAnalog() is specific to the Snow Dots
%   plexonIntegration demo:
%       - it expectes that dotsTheComputablesManager is receiving regular
%       topsDataLog updates, so that the server-side data log matches the
%       client-side data log (i.e., it knows what's going on on the client
%       side).
%       - it expects that the Snow Dots client is sending strobbed words to
%       Plexon's digital inputs, and that each strobbed word has a matching
%       entry in topsDataLog.  It uses these to fit parameters for mapping
%       (a linear transform) Plexon timestamps to client-side topsDataLog
%       timestamps.
%       - it plots the residuals from the clock transformation, showing the
%       relative jitter between strobed word timestamps.  By definition
%       these should cluster near zero.  They should have small variance.
%       - it expects that the Snow Dots client is sending "TTL signals" to
%       the Plexon analog inputs.  The signals should each start with a
%       single rising edge.
%       - it reads data from Plexon's analog inputs, expecting to find
%       signals with a leading, rising edge.  It finds the Plexon time of
%       the rising edge, and maps the time to the client-side topsDataLog
%       frame of reference, using the parameters obtained from strobed
%       words.
%       - it plots the difference between the mapped rising edge time and
%       the clien-side topsDataLog timestamp from the TTL signal output.
%       These should resemble the strobed word residual jitter, but they
%       may be offset and undersampled by the Plexoin analog inputs.
%       .
%
%   If provided, the option argument may be the string "init" or "close".
%   These allow showPlexonAnalog() to initialize and close the Plexon
%   Matlab client SDK, and to set up figures, etc. for plotting.  It
%   expects to be invoked with an "init" and a "close" once each per
%   session.  It expects to be invoked without any option periodically.
%
%   showPlexonAnalog() always returns a result, but it might be empty or
%   meaningless.
%
%   showPlexonAnalog() stores some initialization data, including a Plexon
%   server reference and Matalb graphics handles, in the persistent
%   varible PIDATA.
%

persistent PIDATA
result = [];

if nargin < 1
    option = '';
else
    disp(sprintf('%s: %s', mfilename, option))
end

try
    switch option
        case 'init'
            h.figure = figure(12);
            clf(h.figure);
            
            % raw strobe events and analog signals
            h.strobeRawAxes = subplot(4,1,1, ...
                'Parent', h.figure, ...
                'XTickLabel', []);
            ylabel(h.strobeRawAxes, 'strobed word')
            h.signalRawAxes = subplot(4,1,2, ...
                'Parent', h.figure);
            ylabel(h.signalRawAxes, 'analog signal')
            xlabel(h.signalRawAxes, 'Plexon time(s)')
            
            % strobe times as seen from Snow Dots vs. Plexon clock
            h.strobeAlignAxes = subplot(2,2,3, ...
                'Parent', h.figure);
            title(h.strobeAlignAxes, 'Strobe: plexon timestamps fit to client timestamps')
            xlabel(h.strobeAlignAxes, 'client time(s)')
            ylabel(h.strobeAlignAxes, 'client-minus-Plexon residual(s)')
            h.strobeAlignLine = line(0, 0, ...
                'Parent', h.strobeAlignAxes, ...
                'Color', [1 0 0], ...
                'LineStyle', 'none', ...
                'Marker', '.');
            
            % analog signals aligned using the Snow Dots clock
            h.signalAlignAxes = subplot(2,2,4, ...
                'Parent', h.figure);
            title(h.signalAlignAxes, 'Signal: client timestamp vs. Plexon rising edge time')
            xlabel(h.signalAlignAxes, 'client time(s)')
            h.signalAlignLine = line(0, 0, ...
                'Parent', h.signalAlignAxes, ...
                'Color', [1 0 0], ...
                'LineStyle', 'none', ...
                'Marker', '.');
            
            % various state data
            plexonServer = PL_InitClient;
            PIDATA = struct( ...
                'plexonServer', plexonServer, ...
                'plexonParams', PL_GetPars(plexonServer), ...
                'strobeChannel', 257, ...
                'strobeClientGroup', 'dOut strobed word', ...
                'strobePlexonGroup', 'plexon strobed word', ...
                'signalClientGroup', 'dOut TTL signal', ...
                'signalPlexonGroup', 'plexon TTL signal', ...
                'plexonToClient', [], ...
                'signalFirstChannel', 14, ...
                'edgeThreshold', 1000, ...
                'plotHandles', h);
            %disp(PIDATA)
            
        case 'close'
            PL_Close(PIDATA.plexonServer);
            PIDATA = [];
    end
    
    if ~isempty(PIDATA)
        log = topsDataLog.theDataLog;
        h = PIDATA.plotHandles;
        
        % let axes adjust to new data automatically,
        %   then lock them together, below
        set([h.strobeRawAxes, h.signalRawAxes], 'XLimMode', 'auto');
        set([h.strobeAlignAxes, h.signalAlignAxes], 'YLimMode', 'auto');
        
        % raw strobe events
        [nStamps, timestamps] = PL_GetTS(PIDATA.plexonServer);
        %disp(sprintf('  %d timestamps', nStamps))
        if nStamps  > 0
            strobes = find(timestamps(:,2) == PIDATA.strobeChannel);
            for ii = strobes
                time = timestamps(ii, 4);
                word = timestamps(ii, 3);
                line(time, word, ...
                    'Parent', h.strobeRawAxes, ...
                    'LineStyle', 'none', ...
                    'Marker', '*');
                topsDataLog.logDataInGroup( ...
                    word, PIDATA.strobePlexonGroup, time);
            end
        end
        
        % raw analog signals
        [nPoints, tZero, values] = PL_GetAD(PIDATA.plexonServer);
        nChans = size(values, 2);
        %disp(sprintf('  %d points on %d channels at %.4fs', ...
        %    nPoints, nChans, tZero))
        if nPoints > 0
            params = PIDATA.plexonParams;
            firstChan = PIDATA.signalFirstChannel;
            fChan = params(firstChan + (1:nChans) - 1);
            points = (1:nPoints)-1;
            for ii = 1:nChans
                tChan = tZero + points./fChan(ii);
                vChan = values(:,ii);
                ifFinite = isfinite(vChan);
                col = dec2bin(ii,3)=='1';
                line(tChan(ifFinite), vChan(ifFinite), ...
                    'Parent', h.signalRawAxes, ...
                    'LineStyle', 'none', ...
                    'Marker', '.', ...
                    'Color', col);
                
                % any rising edge?
                %   rise must be at least so large
                risingEdges = find(ifFinite(2:end) ...
                    & (diff(vChan) > PIDATA.edgeThreshold));
                line(tChan(risingEdges), vChan(risingEdges), ...
                    'Parent', h.signalRawAxes, ...
                    'LineStyle', 'none', ...
                    'Marker', '*', ...
                    'Color', 1-col);
                for jj = 1:length(risingEdges)
                    time =  tChan(risingEdges(jj));
                    topsDataLog.logDataInGroup( ...
                        ii, PIDATA.signalPlexonGroup, time);
                end
            end
        end
        
        % lock Plexon time axes together
        strobeLim = get(h.strobeRawAxes, 'XLim');
        signalLim = get(h.signalRawAxes, 'XLim');
        grandLim = [min([strobeLim signalLim]), ...
            max([strobeLim signalLim])];
        set([h.strobeRawAxes, h.signalRawAxes], 'XLim', grandLim);
        
        % strobe timestamp residuals
        if log.containsGroup(PIDATA.strobePlexonGroup) ...
                && log.containsGroup(PIDATA.strobeClientGroup)
            plexonTimes = log.getAllMnemonicsFromGroup( ...
                PIDATA.strobePlexonGroup);
            clientTimes = log.getAllMnemonicsFromGroup( ...
                PIDATA.strobeClientGroup);
            if numel(plexonTimes) >= 2 ...
                    && numel(plexonTimes) == numel(clientTimes)
                tPlexon = [plexonTimes{:}];
                tClient = [clientTimes{:}];
                PIDATA.plexonToClient = ...
                    estimateClockDriftAndOffset(tPlexon, tClient);
                tClientFit = applyClockDriftAndOffset( ...
                    tPlexon, PIDATA.plexonToClient);
                tResidual = tClient - tClientFit;
                set(h.strobeAlignLine, ...
                    'XData', tClient, ...
                    'YData', tResidual);
            end
        end
        
        % signal timestamp-rising edge residuals
        if ~isempty(PIDATA.plexonToClient) ...
                && log.containsGroup(PIDATA.signalPlexonGroup) ...
                && log.containsGroup(PIDATA.signalClientGroup)
            plexonTimes = log.getAllMnemonicsFromGroup( ...
                PIDATA.signalPlexonGroup);
            clientTimes = log.getAllMnemonicsFromGroup( ...
                PIDATA.signalClientGroup);
            if numel(plexonTimes) >= 1 ...
                    && numel(plexonTimes) == numel(clientTimes)
                tPlexon = [plexonTimes{:}];
                tClient = [clientTimes{:}];
                tClientFit = applyClockDriftAndOffset( ...
                    tPlexon, PIDATA.plexonToClient);
                tResidual = tClient - tClientFit;
                set(h.signalAlignLine, ...
                    'XData', tClient, ...
                    'YData', tResidual);
            end
        end
        
        % lock time residual time axes together
        strobeLim = get(h.strobeAlignAxes, 'YLim');
        signalLim = get(h.signalAlignAxes, 'YLim');
        grandLim = [min([strobeLim signalLim]), ...
            max([strobeLim signalLim])];
        set([h.strobeAlignAxes, h.signalAlignAxes], 'YLim', grandLim);
        
        drawnow
    end
    
catch err
    result = err.message;
    disp(err)
end